Source of sdbl.h


/*  sdbl.h  by K.Tsuru  */
// file ID = 30 DRADIX
#ifndef S_DOUBLE_H
#define S_DOUBLE_H
/***************************************************************************
SDouble class
Multi precision real number class in radix DRADIX = 10000

The inner expression.
<<floating point mode(standard form)>>
[+|-]f[0].f[1]...f[ef]...f[s-1]*DRADIX^rdxExp (f[0]=0, f[1]>0)
ef is effective figures. Last h figures are assigned to the hidden precision.
Exponent "rdxExp" has value between INT_MIN (-2147483647 - 1) and INT_MAX 2147483647 on GCC.
<<fixed point mode>>
[+|-]f[0].f[1]...f[ef]...f[s-1]*DRADIX^fixedExp (f[0]=f[1]=...=f[aHead-1]=0)
****************************************************************************/
SDouble SFToSD(const SFraction& x);//SFraction-->SDouble 300
/*------------------------------------------------------------
Radix conversion double x to array res[] in BRADIX.
------------------------------------------------------------*/
int doubleToBinDec(double x, FigBlock& res, uint *size);

class SDouble : public SNumber, public virtual SCalcInfo {
    friend class SDecimal;  // To enable Reform accessible.
    friend class SLong;     // 'figure' accessible in SL::SetSDouble()
// Fixed value of exponent in the fixed point mode.
    static int fixedExp;
// If lower seven bits of fixedPointMode is non-zero, it is fiexed point mode.
// fixedPointMode = depth of nest
    static int fixedPointMode;
// automatic radix conversion BIN_DEC to REAL cannot use.
    SDouble& operator=(const SDecimal& sd);     // has no body
    enum { FIXED_PT_MASK = 0x7f };
    int  rdxExp;        // exponential value in radix = DRADIX
protected:
    void CopyDouble(const SDouble& r, int cs);  // 301
// Functions which set several kinds of value

// By string in type "[+|-]ddd.ddd...[e[+|-]dd]"
    void SetString(const char* s);      // 302
// Substitution double includes integer type variables, long etc.
// If type == BIN_DEC and cyclic decimal, syntax error will occur.
    void SetDouble(double d);   // 303
public:
    void SetLongDouble(long double d);
protected:
    void SetSLong(const SLong& sl); // 304 SLong
// Return an ASCI code of (n+1)-th digit. Uses in function Put().
    int  CharNthDec(long n) const;      // 306
// Make a zero with type = tp.
// Uses in return statement in SDouble.
    friend SDouble SDZero(const SDouble& d);    // 307

// Dec 25, 2000 Move into protected part.
// Temporary maximum size changed by SetEffFig() function.
// Need for the functions using Newton method.
    uint CurrentMaxSize() const {
        return SNEffFig( Type() ) + Hidden() + 1u;
    }
public:
// May 16, 2001 Move into public part from protected one.
// Make a SLong value r = (f[h],...f[t]).
    void ConvertToSLong(uint t, uint h, SLong& r) const;    // 305
/************************************************************************
Additional values for pushCD.
SD_PUSH bit is one, when a SDouble object change the value cutDownTable[].
pushCD = 0...4 are registered in SNumber class, then greater than 4.
**************************************************************************/
    enum { CH_FE = 8, SD_PUSH = 16, TEMP_FREE = 32, REDUCE_MAX_SIZE_DISABLE = 64 };
// set zero(sign = ZERO, all the elements of figure[] are seted by zero)
    void SetZero();     // inline function
/*******************
 Set a short value.
********************/
    void SetInt(int v);         // 308
// rdexp is exponential value in SDouble class
// In SDecimal class rdxExp always set zero for any value of rdxexp.
    void SetSmall(signed char v, int rdexp = 1);
// Effective figures. In BRADIX return a value diffrent from SNManager::EffFigures().
    uint EffFig() const { return SNEffFig( Type() );    }
    uint Size() const {
        return min(SNEffFig( Type() ) + Hidden() + 1u, figure.size());
    }
    uint MinSize() const {
        return (Type() == REAL) ? minArraySize : SNMaxSize(BIN_DEC);
    }
/**********************************************************************
return the exponetial part
[Notice] raw value in fiexed point mode.
The value in the "standard" form i.e. [+|-]f[0].f[1]... * R^e (f[1] >0)
can be obtained by NetRdxExp().
***********************************************************************/
    int  RdxExp() const { return rdxExp; }
    int  NetRdxExp() const {
        return Sign(30) ? rdxExp - (int)aTail +1 : 0;
    }
/******************************************
Set an exponential value in radix = DRADIX
When e = 0, you get the mantissa part.
*******************************************/
    void SetRdxExp(long e) {
        if( (Type() == BIN_DEC) && e ) SetError(SYNTAX_ERR,"SD", 31);
        if(labs(e) > DRADIX_EXP_MAX ) SetError(TOO_LARGE_EXP,"SD", 31);
        rdxExp = (int)e;
    }
/*********************************************************
Return the sign
In the fixed point mode if non-zero elements are out of
maxmum size(CurrentmaxSize()), return zero.
**********************************************************/
    int  Sign(int id) const;    // inline function.
    int  Sign() const { return Sign(32); }
/******************************************************
Return raw sign of SNumber class.
If sign == UNDECIDED, the program abnormally terminate.
*******************************************************/
    int  SNSign() const {   SignCheck(33);  return (int)RawSign();  }
/***************************************************************
Return the decimal exponent in standard form.
                 t
 0.|0000|0000|0aaa|bbbb|...|zzzz| * DRADIX^e
 --> 0.0aaabbbb...zzzz * 10000^(e-t+1)
 --> 0.aaabbbb...zzzz * 10^p   p = 4*(e-t+1) -4 +3
****************************************************************/
    long DExp() const {
        if(Sign(34) == 0) return 0;
        return (long)DFIGURES*(NetRdxExp() - 1L) + iFigures(figure(aTail));
    }
//Return decimal figures
    long DFigures() const;      // 309
//Composed by hidden figures only or not.It is available in fixed point mode only.
    bool IsHidden() const;      // inline function defined bellow
//*this is much less than "a" ?
    bool IsMLT(const SDouble& a) const; // 310
//Return the position of first non-zero element.
//In the standard form aTail = 1 or 0(if sign==0).
//In fixed point mode it is maybe greater than one.
    uint First() const { return aTail < CurrentMaxSize() ? aTail : 0;   }
private:
    uint GetLastPosition() const;       // 311
public:
//Return the position of last non-zero element.
//Last() < CurrentMaxSize()
    uint Last() const {
        if(aHead < CurrentMaxSize()) return aHead;
        return GetLastPosition();   //0.abcd 0000 ..... 0000 1234, return 1
    }
// Head() and Tail() return raw value of SNumber class.
// Please use First()/Last() if possible. ver. 2.17
    uint Head() const {  return aHead;  }
    uint Tail() const {  return aTail;  }
public:
/**************** For making mathematical functions ************************/
/*----------------------------------------------------------------------
Take into the fixed point mode setting the fixed exponent value by "fe".
Use in the calculation of series summantion.
In principle do not call mathmatical functions in the fixed point mode.
The general term of the floating and fixed point modes is named as the
"operation mode".
-----------------------------------------------------------------------*/
    void FixedPoint(int fe);   // 312   DRADIX
/*--------------------
Cancel fixed point mode.
---------------------*/
    //Use with FixedPoint().
    void PointFree();       // 313  DARDIX
    //fixed exponent value
    int  FixedExp() const { return (Type() == BIN_DEC) ? 0 : fixedExp;  }
    //It is in the fixed point mode or not. Return the depth of nest.
    int  IsPointFixed() const {
        return (Type() == BIN_DEC) ? 1 : (fixedPointMode & FIXED_PT_MASK);
    }
    //Take into the fixed point mode for a moment. See Log(x) for usage.
    void TempPointFree(){
        PushCD(PushedCD() | TEMP_FREE);
        fixedPointMode *= (FIXED_PT_MASK+1);
    }
    //Restore the operation mode in one step. Use with TempPointFree().
    void PointModePop() {
        if(!(PushedCD() & TEMP_FREE) ) SetError(SYNTAX_ERR,"PointModePop", 31);
        PushCD(PushedCD() & ~TEMP_FREE);
        fixedPointMode = fixedPointMode/(FIXED_PT_MASK+1) + (fixedPointMode & FIXED_PT_MASK);
    }
/*---------------------------------------------------------------------------
Reform into a representation which is sutable for the present operation mode.
Adjust the size of figure[] i.e. if cutDown != DISABLE it cuts lower needless
zeros.
If the top figure(figure[0]) becomes non-zero by a carry, it shifts to lower
by a figure and makes aTail=1.
Excute after directly changing the elements of figure[] or restoring the
operation mode.
Give to "id" an error ID number of the place in the program.
------------------------------------------------------------------*/
    void Reform(int id);                    // 314
    //Even if in the fixed point mode it changes into a standard form. DRADIX only
    void StdReform(int id);                 // 315
    //Set a long real number by "FigBlock",sign and rdxExp.
    void SetFigBlock(const FigBlock& a, int sgn, int exp){
        SNumber::SetFigBlock(a, sgn);
        SetRdxExp(exp);
        Reform(35);
    }
    //Normalize() version
    void SetFigBlockNorm(const FigBlock& a, int sgn, int exp){
        SNumber::SetFigBlockNorm(a, sgn);
        SetRdxExp(exp);
        Reform(36);
    }
    //If size is greater than "minArraySize" it takes into fixed point mode.
    void FigureAlloc(uint size, int copy){
        SNumber::FigureAlloc(size, copy);
        if(copy <= 0) rdxExp = 0;
    }
/*-----------------------------------------------------------------------------
Try to convert SDouble "d" using long "L" into a form
  d = L * 10^e    ( 0 =< |L| <= ULONG_MAX/DRADIX).
When it is successed, return one. The operation can be reduced into one between
SDouble and short integer.
e.g. d=1.5, x/d ---> (x/15000)*DRADIX.
The multiplication with DRADIX can be done by MultPowRdx(1).
-------------------------------------------------------------------------------*/
    bool ConvTolongExp(long* L, long* e) const;     // 316
/********************************************************************************/

//constructors
    SDouble():SNumber(REAL), rdxExp(0){}
    SDouble(double d):SNumber(REAL), rdxExp(0){
        SetDouble(d);
    }
//  SDouble(long double ld):SNumber(REAL), rdxExp(0){ SetLongDouble(ld); }
    SDouble(const char* s):SNumber(REAL), rdxExp(0){
        SetString(s);
    }
//This is necessary to convert SD*SL into SD*SD.
    SDouble(const SLong& sl):SNumber(REAL){ SetSLong(sl); }
/***************************************************************************
tp=REAL or BIN_DEC
vsz=size of figure[], See FigureAlloc()
If vsz > 0, it is initialized by zero.
If vsz > minArraySize, it takes into fixed point mode and executes
CutDown(DISABLE).
****************************************************************************/
    SDouble(NumberType tp, uint vsz);   //33
// copy constructor
    SDouble(const SDouble& a):SNumber(a), rdxExp(a.rdxExp){}
/*************************************************
Constructor for SFraction --> SDouble.
It will be natural to convert SD*SF into SD*SD.
**************************************************/
    SDouble(const SFraction& sf);       // 34
#if 0
//no use and has no body
    SDouble(const SDecimal& sd);
#endif
//destructor : PointFree()
    ~SDouble();                             // 35
//"SLong = SDouble" has been defined in SLong class.
    SDouble& operator=(const char *s) { SetString(s);  return *this; }
    SDouble& operator=(double d){ SetDouble(d);  return *this; }
//  SDouble& operator=(ldouble ld){ SetLongDouble(ld);  return *this; }  // double version
    SDouble& operator=(const SLong& sl){
        SetSLong(sl);   return *this;
    }
    SDouble& operator=(const SDouble& sd){
        if(this != &sd) CopyDouble(sd, SUBS);
        return *this;
    }
//Convert SFracrion into SDouble.
    SDouble& operator=(const SFraction& sf);    // 36
//Compare absolute values.
    friend int DDCompare(const SDouble& m, const SDouble& n);   // 320
//Auxiliary functions for four operations.
//It cannot use as DsDiv(1.0, 3); for 1/3 = 0.33333... Use a cast as z = DsDiv(SDouble(1.0), 3); since version 2.20
    friend SDouble DDAdd(const SDouble& m, const SDouble& n); // 321
    friend SDouble DDSub(const SDouble& m, const SDouble& n); // 322
    friend SDouble DsMult(const SDouble& m, ulong n);             // 323
    friend SDouble DDMult(const SDouble& m, const SDouble& n);// 324
    friend SDouble DsDiv(const SDouble& m, ulong n);    // 325
//Get the reciprocal of "x" by the Newton's method. It is used in y/x.
    friend SDouble DReciprocal(const SDouble& x);                 // 326
    friend SDouble DDiv2(const SDouble& n);     //diveide "n" by two 327
//four operators
    friend SDouble operator+(const SDouble& m, const SDouble& n);//330
    friend SDouble operator-(const SDouble& m, const SDouble& n);// 331
    friend SDouble operator*(const SDouble& m, const SDouble& n);//324
//Does not provide other operators "SD*int", "SD*long", etc.
    friend SDouble operator*(const SDouble& m, double d);            // 332
    friend SDouble operator*(double d, const SDouble& m){ return m*d;   }
    friend SDouble operator/(const SDouble& m, const SDouble& n);// 333
    friend SDouble operator/(const SDouble& m, double d);            // 334
/*
For two operators "+=" and "-=" when 'n' is equal to zero,
it may be better that it returns "*this" only.
But it is rare case and the copy will cost little time.
*/
    SDouble& operator+=(const SDouble& n) { return (*this = *this + n); }
    SDouble& operator-=(const SDouble& n) { return (*this = *this - n); }
    SDouble& operator*=(const SDouble& n) { return (*this = *this * n); }
    SDouble& operator*=(double d) { return (*this = *this * d); }
    SDouble& operator/=(const SDouble& n) { return (*this = *this / n); }
    SDouble& operator/=(double d) { return (*this = *this / d); }
/**********************************************************************
[Caution]
Please use cast the operations between diffrent type variavles.
SD=0.1*SL yields SD=0.0. --> SD=0.1*(SDouble)SL;
When SL=10, SD=0.1+SL yields SD=10.0 not 10.1. --> SD=0.1+(SDouble)SL;
These are expectative if lhs is SL.
***********************************************************************/
//sign opertors
    SDouble operator+() const { return *this; }
    SDouble operator-() const;  // 335
//multiplied by radix^n
    SDouble& MultPowRdx(int n);     // 340
//multiplied by 10^p, DRADIX only
    SDouble& MultPow10(long p);     // 341
//round off to the last decimal place of effective figures
    enum { STD_REF = 1, CUR_MAX = 2 };
    SDouble& Round(int std = STD_REF);  // 342
/**********************************************************************
Is possible to round off or not.
If the number of continuous (radix-1) or zero from the last decimal place
of effective figures is greater than or equal "n" return 1, else 0.
Use to round off
0.1233 9999 9999 9999 ..... 9999 9875 4578
or
0.1234 0000 0000 0000 ..... 0000 0000 0001
into 0.1234.
Does not include Reform(), then call it if necessary.
[usage]
if(r.IsPossibleToRound(3)) r.Round(0);
************************************************************************/
public:
    bool IsPossibleToRound(uint n) const;   // 343
/****************************************************************************
output functions

0.abcde fghij kl ......... wxyz : perLine
<-fU=5->..................

Put() : not carridge return, Puts() : feed new line
1) If fig(Unit) == 0, it continuously outputs with neither delimiter nor carridge return.
2) putFig : the number of decimal figures.  default : putFig = 0 all figures
3) perLine = the number of blocks which contains "fig" decimal per line
default : perLine = 0 automatically set perLine = displayWidth/(fig+1)
4) lineFormat : NO_CR, CRLF, CONTINUE and END_CR have same meaning as SLong class's ones.
ROUND : do Round() or not.
INT_PUT : It outputs a decimal figure as integer part.
SHOW_FIG: It shows the number of decimal figures at end of line and adds carridge return.
5) Insert delimiter given by "delmt" at intervals of "fig(Unit)" characters.
If delmt == 0, no delimiter is inserted.
It returns the total number of output bytes, including '\n', '.' and delimiter,
as "printf()".
--------------------------------------------------------------------------------------------
void SetFormat(long f, long pf = 0, int pL = 0, int lF = CRLF|ROUND|INT_PUT, int dm = ' ');
******************************************************************************/
    static long figUnit;
    static long putFig;
    static int  perLine;
    static int  lineFormat;
    static int  delmt;
    static long ioCount;    // the number of input or output characters
    static long crCount; // added version 2.30
    enum { NO_CR = 0, CRLF = 1, CONTINUE = 2, END_CR =4 , ROUND = 8,
        INT_PUT = 16, SHOW_FIG = 32 };
    static void SetFormat(long fU, long pF = 0, int pL = 0, int lF = ROUND|INT_PUT|CRLF, int dm = ' '); // 359
    static void DefaultFormat(); // It restores the default values. 360
// They return the number of output charakuters. Please use these to know it in stead of "cout << sd;".
    long Put(long fU=10, long pf = 0, int pL=0, int lF = CRLF|ROUND|INT_PUT, int dm=' ') const; // 344
        // If fU<0 parameters by SetFormat() are used example Put(-1);
    long Puts(long fU=10, long pr = 0, int pL=0, int lF = CRLF|ROUND|INT_PUT, int dm=' ') const;   // x.Puts()
//It outputs the inner expression.
    long RawPut(int crlf = 0) const;                 // 355
    long RawPuts() const {  return RawPut(1);   }
    long CRCount(bool reset=true) const { long temp = crCount;  if(reset) crCount = 0; return temp; }
    friend ostream& operator<<(ostream& os, const SDouble& sd); // added since version 2.21 // 356

//It makes an approximated value taking first "fig" figures.
    SDouble TakeOutFigures(uint fig) const;     // 356
// |*this| == 1 ?
    int IsOne() const {
        if(aHead != aTail) return 0;
        if(NetRdxExp() != 1) return 0;
        return figure(aHead) == 1u ? Sign() : 0;
    }
/*******************************************************************************
It examines whether it is proper to grade up the degree of precision by "upPrec"
figures.When the size of figure[] exceed the power of two only a few, the
efficiency of memory is bad and the merit of FFT multiplication is killed.
It is used in Exp(x), etc.
1.When preferSpeed is ON
 If MaxSize()+upPrec >=2^n, return 0 else 2^n - MaxSize().
 In this case the precision will be lost by few percents.
2.When preferSpeed is OFF(addvance precison,default)return "upPrec".
*******************************************************************/
    uint ProperUpPrec(uint upPrec) const;       // 357
}; 

/**************** end of SDouble class ******************/

//It moves "fixedPointMode" value to upper bits.

inline int SDouble::Sign(int id) const {
    int snsign = SNumber::Sign(id);
    if(!snsign || !IsPointFixed() || (Type() == BIN_DEC)) return snsign;
    int e = (int)aTail - rdxExp + fixedExp;
    return (e < (int)CurrentMaxSize()) ? snsign : 0;
}
inline void SDouble::SetZero(){
    rdxExp = 0;     SNumber::SetZero();
}
inline long SDouble::Puts(long fig, long pr, int perLine, int lF, int delmt) const {
    long l = Put(fig, pr, perLine, lF, delmt);
    FPutc('\n');
    return l+1;
}
inline bool SDouble::IsHidden() const {
    if(SNumber::Sign(37) == 0) return true;
    return aTail > EffFig() ? true : false;
}
// SDouble contants added since version 2.30
static const SDouble ONE(1.0);  // one
static const SDouble ZERO(0.0); // zero

/*  defined in template in "mathtp.h"
inline istream& operator>>(istream& s, SDouble& a) // added since version 2.21
{
    string buf;
    std::getline(s, buf);   // usage : cin >> a; ( s = cin )
    a = buf.data();
    return s;
}
*/
#endif  // S_DOUBLE_H
sdbl.h : last modifiled at 2017/10/18 15:19:28(19,674 bytes)
created at 2016/04/11 11:18:58
The creation time of this html file is 2017/10/23 10:29:22 (Mon Oct 23 10:29:22 2017).